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Abstract 

Consider the Laplacian in a bounded domain in M. d with gen- 
eral (mixed) homogeneous boundary conditions. We prove that its 
eigenfunctions are 'quasi-orthogonal' on the boundary with respect 
to a certain norm. Boundary orthogonality is proved asymptotically 
within a narrow eigenvalue window of width o(E 1 ^ 2 ) centered about 
E, as E — > oo. For the special case of Dirichlet boundary conditions, 
the normal-derivative functions are quasi-orthogonal on the boundary 
with respect to the geometric weight function r • n. The result is in- 
dependent of any quantum ergodicity assumptions and hence of the 
nature of the domain's geodesic flow; however if this is ergodic then 
heuristic semiclassical results suggest an improved asymptotic esti- 
mate. Boundary quasi-orthogonality is the key to a highly efficient 
'scaling method' for numerical solution of the Laplace eigenproblem 
at large eigenvalue. One of the main results of this paper is then to 
place this method on a more rigorous footing. 



1 Introduction and main results 

Let Q C M. d be an open bounded Euclidean domain of dimension d > 2, 
with boundary T = dfl parametrized by the (d— 1) dimensional coordinate 
s G T, and outwards unit normal vector n(s). For instance we envisage the 
boundary being a finite union of smooth surfaces with angles at junctions 
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bounded away from zero and 2ir; in general our domain can be Lipschitz. 
The set of eigenfunctions {0j(r)}, i — 1, 2 . . . oo of the Laplace operator A 
are defined by 

-Afa = Ei<j>i, (1) 

everywhere in Q, and by Dirichlet boundary conditions on some subset of 
the boundary, 

<f>i(s) = for s G T D , (2) 
and homogeneous boundary conditions on the remainder, 

7(8)0,(8) + d n( f )i (s) = forser\r D , (3) 

where <9 n 0j(s) = n(s) ■ V0i(s) is the normal derivative at the boundary. 
We assume the subset T D C T is composed of a finite number of compact 
pieces, and that on the remainder 7(s) is bounded. The boundary conditions 
are such that A is self-adjoint, and we also assume they are such that a 
complete set of eigenfunctions exists. The corresponding eigenvalues are 
ordered E\ < E 2 < • • • , and eigenfunctions are orthonormal on the domain, 
(<f>i,(f>j)n = Sij. In keeping with common quantum-mechanical terminology 
we will say that level % has energy Ei. 

Let q be the following symmetric bilinear form depending on an energy 
parameter £, 

q(u,v;£) := j> — - — (u n v + v n u) + r n ^uv - Vu ■ Vv^j + u r v n + v r u n dA. 

(4) 

Its particular form has properties that will become apparent shortly. Here 
u and v are functions defined with their derivatives on fi, the closure of f2. 
We use the abbreviations u n = d n u, u r = d T u = r • Vw, and r n = r ■ n. The 
surface element on V is dA. Note that q is sensitive only to boundary values 
and first derivatives of u and v. 

Let Q be the semi-infinite symmetric matrix defined using the eigenfunc- 
tions of domain Q by 

Qij ■= q(<i>i,<i>j;£ij), i,j>l (5) 

where = Ei + Ej is the sum of the energy eigenvalues. 
Our main result is 
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Theorem 1 (Quasi-orthogonality) 

onal elements 

Qii = 2£/j, 

and off- diagonal elements satisfying 
\Qij\ < C(Ei — Ej) 2 , 



The matrix Q defined above has diag- 
for all i > 1, (6) 

for alii ^j, i,j>l, (7) 



where the constant C depends only on the shape of the domain. 

Note that by (6) and (7) the diagonal elements grow linearly with energy 
whereas off-diagonal elements are bounded by a fixed quadratic function of 
the energy difference. 

Bearing in mind that eigenfunction derivatives grow like ||V0i||x,2(n) = 

1 /2 

0(E i ), we might naively expect from the size of terms in (4) that all off- 
diagonal elements of Qij near the diagonal (Ei Ej) grow linearly with 
Ei, as Ei (and therefore Ej) tends to infinity. Theorem 1 tells us that this 
is not so, and that the choice of q carries with it non-trivial cancellation. 
This particular property of q motivates our choice of the bilinear form. We 
conjecture that this property is unique to the form q given by (4), up to 
a choice of origin. It should be emphasized that this result holds for all 
elements not just asymptotically. We immediately have the following 

Corollary 1.1 (Asymptotic orthogonality) Given Q as defined above, 
let Q be the sub-matrix of Q corresponding to all levels i,j in a local energy 
window Ei,Ej e [E — €q, E + eo] where e = 0(E lS ). Then for all (3 < 1/2, 

1. 

as e^oo, (8) 

where I is the identity matrix. Thus, close to the diagonal, boundary 
orthogonality is asymptotically exact. 

2. Off-diagonal elements ofQ/2E are bounded in absolute value by cpE 213-1 

The sub-matrix Q is positioned in Q as indicated in Fig. la. Note that 
by Weyl's law for the asymptotics of the level counting function [16], #{i : 
Ei < E} ~ c d vol(Q)E d / 2 , it follows that the limit E — > oo is equivalent to 
i,j — > oo. It also follows that for all d > 2, a growing number of levels 

~ Q vol(Q)£ , ^ _1 )/ 2 falls within an energy window e ~ E 1 / 2 . This means 



Figure 1: a) Illustrates the extraction of Q from the full matrix Q. b) 
Illustrates (for d—2) the geometry of the domain, and its escribed circle. 

that, since (3 can take any value strictly below 1/2 that the order N of the 
matrix Q which tends to the identity can be allowed to grow without limit 
like N = o(N ). 

We now focus on the special case of Dirichlet boundary conditions, 

<MH = o, (9) 

corresponding to T D = T. In this case a simple calculation shows that the 
quadratic form (4) simplifies such that (5) becomes 

Qij = j>r n d n (t>id n (t)jdA, (10) 

that is, a boundary inner product under a weighting function r n . In this case 
Corollary 1.1 can then be restated as 

Corollary 1.2 (Dirichlet quasi-orthogonality) As E — > oo the eigen- 
function normal derivatives belonging to an energy window of width o(E l t 2 ) 
become asymptotically orthogonal under the boundary weighting function r n . 

This Corollary is a stronger form of a conjecture first made (we believe 1 ) 
by Vergini-Saraceno (see Appendix of [25]). Making use of an identity equat- 
ing (10) with (Ej - f?i)(0i,r • V0j) n for % ± j (e.g. see [4] Eq.(H.25)), they 

1 Although it was known [8] to be exact for the degenerate case Ei = Ej. 
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Figure 2: Density plot of the matrix Qij for the range i,j = 1---N for 
N = 1000. % — j — 1 is in the top left corner. The domain used was a smooth 
deformation of the circle, in d — 2, with Dirichlet boundary conditions; the 
classical flow is almost completely ergodic. By Weyl's Law the energy level 
density is constant, therefore the axes can equally well be thought of as Ei 
and Ej. Eigenfunctions were calculated using the scaling method ([5], Ap- 
pendix B). Dark indicates large magnitude values of Q, white indicates zero. 
The visible white band (very small values) with growing width E 1 ^ 2 around 
the diagonal (dark line) shows the quasi-orthogonality property special to 
the choice of (4). The visibly 'banded' structure is discussed in Section 3.1. 
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suggested that off-diagonal elements vanish linearly as \Ei — Ej\. We now 
know by Theorem 1 that the vanishing is of higher order, namely as (Ei—Ej) 2 . 
Understanding this order of vanishing is important for understanding the ac- 
curacy of the numerical 'scaling method', as we will discuss below. 

For domains whose geodesic flow (equivalently, classical dynamics) is er- 
godic, a semiclassical argument (of the type we will discuss in Section 3.1), 
more powerful than that of Vergini-Saraceno, was developed by Cohen, Heller 
and this author. This predicted the (Ei — Ej) 2 factor appearing in (7), and 
showed excellent agreement with numerical studies [3, 4]. The argument in- 
volved evaluation along classical trajectories of a non-smooth operator which 
'lives' (is a delta-function) on the boundary and possesses a 'strength' D(s) 
(see [4]). For the case D(s) = r n , corresponding to dilation of the billiard, the 
operator can be shown to be an exact second time- derivative of a bounded 
function. As a result its power-spectrum vanishes at zero frequency and this, 
via the semiclassical argument, is seen to be the cause of the (Ei — Ej) 2 factor. 
However this explanation has problems: i) it relies on a semiclassical argu- 
ment not yet proven to hold for all matrix elements (see Section 3.1), ii) the 
semiclassical argument is not known to be valid for such singular operators 
living on the boundary, and iii) only ergodic domains can be handled. 

Corollary 1.2 provides proof of Dirichlet quasi-orthogonality without re- 
sort to any assumptions about the shape or ergodicity of the domain. In effect 
the identity in Lemma 1.1 below bypasses the two classical time-derivatives, 
and the semiclassical argument, by directly accounting for the (Ei — Ej) 2 
factor. Thus quasi-orthogonality is seen to be a result independent of any 
semiclassical or quantum ergodic assumptions. What we will be left with 
is consideration of matrix elements of the operator r 2 , which is a smooth 
and bounded function of space. Bounds on these matrix elements are much 
easier to find than in the case of matrix elements involving boundary infor- 
mation. Lemma 1.1 relies on overlap identities which follow from the Diver- 
gence Theorem alone; these identities are derived in the Appendix using a 
symbolic matrix technique. Theorem 1 will then follow from an elementary 
bound on matrix elements. The possibility of considering boundary quasi- 
orthogonality with general boundary conditions, as we do with Theorem 1, 
is to our knowledge new. 

In terms of applications, Dirichlet quasi-orthogonality (Corollary 1.2) is a 
key component of the 'scaling method' invented by Vergini-Saraceno [25] (see 
[4, 6], and Appendix B of the companion paper [5] for review), a very efficient 
method for the numerical solution of the Dirichlet eigenproblem in strictly 
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star-shaped domains (r n > 0). This method has enabled large-scale studies 
of eigenfunctions at extremely high E in d = 2 [24, 5] and d = 3 [20], where 
it outperforms all other known methods by a factor of typically 10 3 . The role 
played by the matrix Q is as follows: The small size of off-diagonal elements of 
Q near the diagonal means that a set of the domain's Dirichlet eigenfunctions 
approximately diagonalize two certain quadratic forms (see Appendix B of 
[5]). The numerical method involves simultaneous diagonalization of these 
two quadratic forms, in order to compute a large number of approximate 
domain eigenfunctions at once. The error in the approximation, and therefore 
in the functions found by the method, depends on among other things the 
order of vanishing of Q as one approaches the diagonal. Thus our work 
finally places this important method on a more rigorous footing, one which, 
despite its successful use, it has not yet had. In particular its success in 
domains with mixed dynamics (divided phase space), quasi-integrable, or 
integrable dynamics is now no longer mysterious. Moreover the generality of 
(2) and (3) opens up the possibility of extending the scaling method to solve 
eigenproblems with more general boundary conditions. 

As another application, in the Dirichlet case (6) becomes 

j r n {d n <pi) 2 dA = 2E h for all i > 1, (11) 

a result first found by Rellich [21] and independently by Berry- Wilkinson [8]. 
(Other more simple derivations have recently been found [3, 4, 13]). This for- 
mula is useful numerically because it enables Dirichlet eigenfunctions to be 
correctly normalized using boundary information alone. For general bound- 
ary conditions (6) reads 

j{d - 2)0A0* + r n (Erf? - |V0i| 2 ) + 2d r &d n & dA = 2E t , for all i > 1, 

(12) 

a normalization formula involving boundary values and first derivatives alone, 
which we have not found in the literature. By contrast, the general boundary 
condition formula of Boasman for d = 2 (Appendix D of [9]) requires second 
derivatives, a numerically more complicated demand. 

Finally we note some other recent work on eigenfunctions on the bound- 
ary. In the Dirichlet case the asymptotic completeness in C°°(r) of the 
boundary functions has been shown (this holds even when restricted to an 
energy window), and the curvature correction to their local intensity derived 
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and tested numerically [2]. A version of the Quantum Ergodicity Theo- 
rem [23] for boundary functions has been proved in piecewise-smooth do- 
mains [14]. We note that boundary normalization formulae such as (11) 
break down for general Riemannian manifolds with boundary [13], so we 
expect that our results are specific to the Euclidean case. 

The remainder of the paper is structured as follows. In Section 2 we prove 
Theorem 1, and give a recipe for translating the origin in such a way that 
the bound C achieves its minimum value Co- In Section 3 we review some 
semiclassical results on the size of matrix elements of well-behaved operators 
as E — > oo. Applying them to the operator r 2 , for domains with ergodic 
flow we improve somewhat the convergence rate in part 2) of Corollary 1.1. 
However, in contrast to our main results, this improved ergodic statement is 
not proved to hold for every single element i,j. We will also briefly mention 
expected behaviour for domains with integrable and mixed flow. The Ap- 
pendix contains a detailed account of a general symbolic matrix procedure 
used to derive certain identies used in the proof in Section 2. 

Acknowledgements. This work could not have been possible without 
interactions with and important feedback from Percy Deift, Fanghua Lin, Pe- 
ter Sarnak, Maciej Zworski, Kevin Lin, Doron Cohen and Eduardo Vergini. 
The boundary formulae of Appendix A originated in work with Michael Hag- 
gerty and Eric Heller. The author is supported by the Courant Institute at 
New York University. 

2 Proof of Theorem 1 

Theorem 1 hinges on the following identity between elements of the matrix 
Q and matrix elements of the function r 2 over the domain. 

Lemma 1.1 Let Q be the matrix defined by (4) and (5), with Laplace eigen- 
functions defined by (1), (2) and (3). Then, 

Qij = 2E i 5 ij + ^-(0;,r 2 0,,)n, for all i,j > 1, (13) 

where €ij = Ei — Ej is the difference in energies. 

Notice that in (13) the first (resp. second) term contributes only on (resp. 
off) the diagonal. 



8 



Proof. Consider functions u and v, defined with their derivatives on Q, 
and each satisfying the Helmholtz equation, namely, 

— Au = E u u, 

-Av = E v v, (14) 

everywhere in Q. The functions need satisfy no particular boundary condi- 
tions on T. There are two cases which are handled separately: the energies 
E u and E v are either equal or unequal. Examining first the equal-energy case 
E u = E v = E, in the Appendix we derive the identity (see (36)), 

(u,v) n = -^j> —^—( u nV + v n u) + r n (Euv - Vu ■ Vv) + u r v n + v r u n dA. 

(15) 

This identity is a consequence of the Divergence Theorem applied to various 
combinations of the functions u and v and their first derivatives. Comparing 
this with (4) gives 

(u,v) n =^q(u,v;2E). (16) 

Substituting u = <pi and v = <f>j and using the (5) and orthonormality on the 
domain gives, 

Qn = 2E i 6 ij , for Ei = Ej. (17) 

This covers both the diagonal element case i = j and the degenerate case 
where i ^ j but Ei = Ej . 

Now examining the case E u ^ E v , in the Appendix we derive the identity 
(equivalent to (33)), 



e 2 , , . fd-2, , E 



e 



^{u,r v) n = j)^—(u n v + v n u)+y---rj(u n v-v n u) 

+ r n (^ uv — • Vf ^ + u r v n + v r u n dA, (18) 

where the energy difference is e = E u —E v and the sum S = E u +E v . Applying 
self-adjoint boundary conditions (2) and (3) causes the term anti-symmetric 
in u and v to vanish identically on the boundary, leaving 

— (u,r 2 v) n = q(u,v; E u + E v ). (19) 



9 



Substituting u = (pi and v = <f>j and using (5) gives 

Qij = ^-{hytih for£,^ 3 , (20) 

The two cases — equal energy (17) and differing energy (20) — can be summa- 
rized by the one identity (13). □ 

We now see that the particular form (4) of q results from a relatively 
complicated sequence of manipulations detailed in the Appendix. Lemma 1.1 
has reduced the question of quasi-orthogonality to one of the size of the ma- 
trix elements ((pi,r 2 (pj)n. Now the diagonal elements (0j,r 2 0j)n are trivially 
bounded by 



r z 4>idV < R 2 / (pjdV = R z , (21) 
n Jn 

where dV is the volume element, and the maximum radius from the origin 
attained by the domain is 

R = maxlrl. (22) 
rer 

Then by Cauchy-Schwarz we have off-diagonal matrix elements bounded by 
the same constant, 

m,r%)n\ < |(^,r 2 ^)nr /2 |(0i,r 2 ^)nr /2 < R 2 ■ (23) 

Substituting this into (13) completes the proof of Theorem 1, with the con- 
stant being C = R 2 /A. □ 

Clearly the eigenproblem (1), (2) and (3) is translationally invariant, that 
is, {Ei} and {(pi} are independent of the choice of origin for our coordinate 
r. It follows from (6) that the diagonal elements of Q are also translationally 
invariant; however the off-diagonal elements of Q are not. Let us assume 
our task is, given a domain Q of particular shape, to find the boundary form 
given by (4) for which the constant in (7) can be chosen to be smallest. We 
are free to translate the origin to achieve this goal. Clearly the minimum C is 
given by Cq = -Rq/4, where Rq is the radius of the smallest circle (generally, 
ball in IR d ) enclosing Q. This is the escribed circle (or ball); see Fig. lb. The 
optimal choice of origin is the center of this circle, ro- We note that any 
choice for the origin need not fall inside Q, and that none of the results given 
in this paper depend on whether it does or not. (This contrasts the scaling 
method, for which it is believed that the domain must be strictly star-shaped 
with respect to the origin). 
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3 Semiclassical estimates 



The power of Lemma 1.1 is that it has reduced the quasi-orthogonality ques- 
tion to one of the size of matrix elements of a bounded operator r 2 , in such 
a way that a trivial bound on ((f)i,r 2 (f)j)ci is adequate to prove Theorem 1. 
Can we go beyond this bound? Estimating matrix elements semiclassically 
(at large E) has been a major theme in both the physics and mathematics 
communities, and remains an active area of research. We will now draw on 
their results. We emphasize that to prove Theorem 1, semiclassical results 
have not been used. 

Depending on the domain ('billiard') shape, the geodesic flow on Q x S* 6 ^ 1 , 
that is, the classical motion of a trapped point particle undergoing elastic 
reflections from T, falls into three broad categories: ergodic, integrable, and 
mixed [18]. This has consequences for the Laplace eigenf unctions which have 
been a major theme in 'quantum chaos' [16, 22]. 

3.1 Ergodic domains 

In the ergodic case, the Quantum Ergodicity Theorem [23] (QET) has been 
proved, stating that asymptotically (E — > oo) all but a set of measure zero 
of the eigenfunctions become spatially equidistributed across the domain. 
Thus almost all diagonal elements (0j,r 2 0j)Q tend to the classical average 
f n r 2 dV /vol(Q). However this is not of direct use since it is only off-diagonal 
elements of r 2 that play a role in Lemma 1.1. Assuming a statistical model in 
which eigenfunctions behave like uncorrelated random- waves [7] , off-diagonal 
elements are expected to be Gaussian distributed with zero mean; this has 
indeed been verified numerically [1, 19]. 

The remaining issue is their variance. This issue is discussed in much 
more detail in [5] but we provide a summary here. Almost all off-diagonal 
elements have been proven to vanish asymptotically [27] implying that the 
variance must vanish. However, the random-wave prediction for variance 
size is poor [1, 3, 5] since it cannot account for variance changing as a func- 
tion of energy difference, that is, a 'banded structure' to the matrix. The 
banded matrix structure in the case of an ergodic domain is shown in Fig. 2. 
The semiclassical sum rule of Feingold-Peres [12, 11] (FP) has proven much 
more successful in predicting numerically observed variance [12, 11, 1, 19, 3], 
including the banded structure. There is numerical evidence [5] that it is 
asymptotically correct in billiard systems, but that convergence is quite slow. 
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The band profile of the matrix elements of an operator A, that is, variance 
as a function of energy difference, is related to Ca(w), the Fourier trans- 
form of the auto-correlation of the corresponding classical operator under 
the (unit-speed) classical flow. FP states that for operators A whose classi- 
cal counterparts are well-behaved in phase-space, 



var 



(0,, Mj)a] - ^0l E ~ 1 ' 2 for E if *E jt »E^ oo, (24) 



where the wavenumber difference is u>ij = E X J 2 — E 1 ^ 2 . The sense in which 
we use 'variance' here is the variance of a large sample of matrix elements 
which share similar values of lo^ and fall within a classically-small energy 
range. Thus, we are measuring variance within a small 'patch' of the matrix 
A(j)j)n. As with Corollary 1.1, we now restrict ourselves to the energy 
window e = o(E 1 ^ 2 ), within which ofy- — > as E — > oo (corresponding 
classically to taking the zero-frequency limit). Choosing A = r 2 , then (24) 
and Lemma 1.1 gives matrix element size 

4^olf| J E-^M-E*) 2 ' ( 25 ) 

which holds for nearly all off-diagonal elements, and should be compared to 
the hard bound (7). The domain-dependent constant (7 r 2 (0) can be estimated 
numerically via trajectory simulations; we do not know if it is minimized by 
the optimal choice of origin r derived in Section 2. Thus in part 2) of 
Corollary 1.1, we can improve the bound to cg-E 2/3_5//4 . For (3 = (constant 
energy difference), the convergence rate then improves somewhat from E^ 1 
to E~ 5 / 4 . Note that this result is equivalent to that of [3, 4], but because 
r 2 is a smooth rather than singular function the approach presented here 
stands on more solid foundations, as discussed in Section I. However, FP 
is not proven to hold for every single matrix element: there may exist an 
exceptional set (analogous to the diagonal case) which fails to converge as 
above. Theorem 1 remains the only proven bound on Q known to us. 

It is worth pointing out that ergodicity and the FP sum rule, if it were 
valid for singular functions (as it empirically appears to be [3]), would imply 
a weaker form of asymptotic quasi-orthogonality, regardless of the form of q. 
That is, Corollary 1.1 would hold even for generic boundary forms with terms 
similar to those in (4) but lacking their particular cancellations, for almost 
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all off-diagonal elements, independent of Ei — Ej, but with the much slower 
convergence rate ~ cE~ 1 ^. In the case of Dirichlet boundary conditions, 
this would correspond to choosing a generic weight function D(s) in place of 
r n in (10). 

When the ergodic flow furthermore is Anosov (uniformly hyperbolic), a 
stronger version of QET has been conjectured [22], which lifts the require- 
ment that a zero-density set be excluded. This is known as Quantum Unique 
Ergodicity (QUE), and currently remains unproven in Euclidean domains, 
although good numerical evidence exists when d — 2 [5]. Since QUE can 
be extended to off-diagonal elements [28], if proven, QUE would enable the 
improved bounds (25) to be claimed for every single off-diagonal element in 
these classes of domains. 

If the flow contains orbits with zero Lyapunov exponent ('bouncing ball 
modes'), classical autocorrelations generally have a power-law rather than 
exponential tail, thus Ca{w) diverges as uj — > 0, and (25) is undefined. In 
this case FP must be modified, and to our knowledge only the diagonal 
variance has been carefully studied theoretically [11]. 

3.2 Integrable and mixed domains 

Finally we briefly mention some of what is known about the behaviour of 
matrix elements in systems with other categories of classical dynamics, and 
speculate on consequences for the rate of convergence in Corollary 1.1. 

If the classical dynamics is integrable, then eigenfunctions are regular, 
with Wigner functions concentrating around d- dimensional classical invari- 
ant tori in phase space [18, 16]. Each eigenfunction possesses a well-defined 
set of quantum number labels p G Z d ; values of r tend to be uncorrelated as a 
function of level number. Much less has been said about off-diagonal matrix 
elements in this case. In smooth Hamiltonian systems matrix elements be- 
tween levels labelled by p and q are expected to die like exp(— c|p — q|), for a 
generic smooth operator. Thus most off-diagonal elements are exponentially 
small (the so-called 'selection rules' [16, 17]), but a few may be large, that is, 
of the same order as diagonal elements. However, because our domain has a 
hard boundary condition, it is easy to estimate (by analogy with the d — 1 
case) that generic matrix elements die algebraically with |p — q|). Thus we 
still expect most off-diagonal elements to be small. This has been observed 
numerically in billiard systems [19, 17]. Thus, we expect that convergence to 
quasi-orthogonality will be much more rapid, for the majority of matrix ele- 
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ments, than for the ergodic case. In general, it is hard to say more. It would 
be interesting to see if integrable examples can be found where (0j,r 2 0j)n 
approaches R 2 as E — > oo; this would prove that Theorem 1 is sharp, with 



A generic domain has mixed dynamics, and eigenfunctions either occupy 
ergodic components of phase space, or lie on integrable tori. Matrix ele- 
ments involving energy levels lying in different components are therefore ex- 
pected to be very small. This has been verified numerically in some smooth 
Hamiltonian systems [10], but no study in billiards is known. There is much 
opportunity for numerical study in such billiard systems. 

Appendix A: Boundary identities for matrix 
elements 

We present and apply a novel general procedure for deriving identities which 
express certain bilinear forms over a domain, involving solutions to the 
Helmholtz equation purely in terms of bilinear forms on the boundary. Given 
a bounded open domain Q C M. d , d > 2, let the fields u(r), v(r) be regular 
solutions of (14) with energies E u and E v respectively, everywhere in Q. We 
envisage T being the union of a finite number of smooth curves (surfaces), 
however the broader condition that Q be Lipschitz is sufficient. Crucially, we 
do not require any particular boundary conditions to be satisfied by u, v on 
T. We will make use of the Divergence Theorem applied to vector fields which 
are bilinear combinations of u, v and their first derivatives; thus u and v need 
only be sufficiently regular that the Divergence Theorem can be applied to 
such a vector field. These conditions on vector fields are quite broad [15]. 
We believe that the conditions laid out at the beginning of this paper are 
sufficient to ensure that u or v can be set to equal any eigenfunction 4>i and 
still have the Divergence Theorem apply. We will express 



dV being the volume element in R d , in terms of boundary integrals involving 
values and first derivatives of u and v, for the operator choices H — 1 and 
H = r ■ r := r 2 . We will consider both the cases E u ^ E v and E u = E v . 

Consider functions on Q that are built from bilinear forms in u, v and 
their spatial derivatives. We restrict attention to such functions that are 



C = R 2 



/4- 




(26) 
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order 


r o 






go 


uv 






d 1 


uv 


ruv, <-> 


fruv, rurv 


d 2 


uv 




^ ^N^-N^ / S . 1 

rruv, ruvr, rurv, <-> 



Table 1: Lowest order scalar diagrams ('molecules') from which the set of 
functions {qp} was chosen. An order of d m means that m is the highest 
order derivative of u or v. The symbol <-> means a repetition of the previous 
diagram but with -u and v swapped. See text for further explanation. 



order 


r o 




^2 


go 








d l 


iiv, <-> 


™, rm>, <-> 


fruv, 





Table 2: Lowest order vector diagrams ('radicals') from which the set of 
functions {p a } was chosen. See Table 1 caption for explanation. 

either scalar fields or vector fields. The key is to select a set of scalars {qp}, 
(3 = l,2...N q , and a set of vectors {p a }, a = 1,2...N P , such that the 
divergence of each vector field can be written as a (location-independent) 
linear combination of the scalar fields. Applying the Divergence Theorem 
turns these into relations between boundary and volume integrals. The coef- 
ficients form a matrix, which can be handled symbolically to express volume 
integrals as linear combinations of boundary integrals. We note that the 
idea of inverting such a coefficient matrix originated in unpublished work of 
Michael Haggerty (see Appendix H of [4]), however here we extend the idea 
and explore a larger set of vector and scalar terms. 

To build the functions we use the following objects: the fields u, v, their 
gradients Ui := d(u and higher derivative tensors U; L j := didjU, etc, and 
polynomials in the coordinates r\ which transform as scalars, vectors, tensors, 
etc. It is possible to enumerate systematically, using diagrams, the hierarchy 
of such combinations with overall scalar or vector character. Diagrams also 
greatly ease the algebra in taking the divergence of the vector functions (we 
apply the rules by hand, but automating them would not be hard) . 

We have found the following scheme very useful. A combination of overall 
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scalar character is represented by a 'molecule' (or a disconnected collection of 
molecules) built from 'atoms' connected by 'bonds'. Each atom can be either 
r, u, or v. A bond carries summation over a spatial index, for instance i; 
when a bond terminates at u or v this corresponds to a derivative di, whereas 
termination at r corresponds simply to rj. Thus u and t> may carry 0, 1, . . . oo 
bonds whereas each r must carry exactly 1. See Table 1 for the few lowest- 
order scalar functions. For example, fuv represents Yli=i Sj=i r % u %j v y From 
now on summation over repeated indices is implied. 

A combination of overall vector character is similar except it has exactly 
one 'dangling bond', represented by a dot (we could think of it as a 'radical'). 
Taking the divergence of a vector function corresponds to first summing the 
molecules resulting from connecting the dangling bond to each of the atoms 
in turn (including the one from which it originates). The molecules produced 
are not necessarily valid, so the following three 'reaction' rules need to be 
applied repeatedly until each molecule becomes valid: 

1. u — > —E u u. That is, self-bonds to atoms u or v vanish and are replaced 
by the scalar multiplier — E u or — E v respectively. 

2. r — > d. That is, any r bonded to itself vanishes along with the bond, 
and is replaced by the scalar multiplier d. 

3. IT — > ~. That is, any r with two bonds (not covered by the previous 
rule) vanishes, and the resulting dangling ends connect to form a single 
bond. 

It can be easily verified that this procedure is equivalent to taking the di- 
vergence of any given vector field. For example, the divergence of r 2 UiV is 
computed as follows, 

rriiv — > 2rruv + fruv + fruv 

r ^ s 2ruv — E u rruv + rriw, (27) 

thus the answer is 1r\UiV — E u r 2 uv + r 2 UiVi, evident as row a = 7 of the 
matrix equation (28) below. 

Using this technique we have calculated the divergence of the simplest 25 
or so vectors, which results in a similar number of scalars. For the goal at 
hand we have eventually found the following vector subset of size N p = 8 
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and scalar subset of size N q = 8 to be minimally sufficient, 
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This set of divergence relations can be more compactly written, 

N q 

V-Pa = ^Mapqp, a = l,2...N p . 

13=1 

Integrating over Q and applying the Divergence Theorem gives 

N„ 



(28) 



(29) 



fn-p a dA = y^Mafs / qpdV, 
Jr p =1 Jn 



a = 1, 2 . . . 7V p . 



(30) 



We first consider the case E u ^ E v , for which the matrix M is invertible 
(in fact we chose the above set of vectors to make this so). Thus each 
volume integral appearing in (30) can be written purely in terms of boundary 
integrals, 



r Np r 

q a dV = ^(AT 1 ) fn-ppdA 

Jn 0=1 Jr 



a = l,2...N q . (31) 



A symbolic algebra package (Mathematica) was used to find .A/f -1 , which 
we need not write out in full here. Rather we use only select rows of 7W x , 
that is, certain values of a in (31). For example, the a = 1 row of 
is (1/(E U — E v ), — 1/(E U — E v ), 0, 0, 0, 0, 0, 0), immediately gives the simple 
relation 

(u, v) Ua = — — i {un -Wv-vn- V«) dA, (32) 

E u — E v J Ta 
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well-known from Green's Theorem. Choosing the a = 7 row gives, expressed 
in terms of e := E u — E v and £ := E u + E v , the relation 

J r uvdV = — J -^—{u n v + v n u) + I - - -r I (u„f - v„w) 

+ r n ^7^ m ' — ^ M ' Vl^ + M r^ri + ^rW„ cL4- (33) 

In the particular case of Dirichlet boundary conditions on T, using Vw = nu n , 
etc, this relation simplifies to 

/ r 2 uvdV = (fr n u n v n dA, (34) 
./n e ./r 

a result we have not found in the literature. 

We now consider E u = E v = E, in which case M is not invertible (its 
first two rows are identical). However J n q a dV can still be expressed as 
a linear combination of boundary integrals if the standard unit basis vec- 
tor := (0, . . . , 0, 1, 0, . . . , 0) G M. Nq , which contains a 1 only in loca- 
tion a, lies in RowA4. Since in our case Nul.M contains the single vector 
(0, 0, 0, 0, 0, 0, 1, E), this can be done for a G [1, 6]. For each such a, we can 
find the row coefficient vector x( Q ) G M. Np by solving the (consistent and in 
this case underdetermined) linear equation 

M T x (a) = b (a)_ ( 35 ) 

We remind the reader that the coefficient vector spaces between which Ai T 
transforms should not be confused with coordinate vectors in M. d . The so- 
lution set of (35) is readily found with a symbolic algebra package. For the 
case a = 1 the resulting row coefficients give the identity 

f 1 f d — 2 

/ uv dV = — — <£ (u n v + v n u) + r n (Euv — Vw • Vi>) + u r v n + v r u n dA. 

J n 2E J r 2 

(36) 

where we have (arbitrarily) chosen to equalize the first two coefficients in 
order to make the u^v symmetry manifest. This identity involves no deriva- 
tives higher than the first, and to our knowledge is new. A similar identity for 
d — 2 which required second derivatives, and is therefore more cumbersome, 
has already been found [9] (also see [4], Eq.(H.7)). In the particular case of 
Dirichlet boundary conditions (36) directly proves (11). 
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What could be gained by enlarging the function spaces beyond N p = 8 
and N q = 8 ? For both cases E u = E v and E u ^ E v , the only criterion 
for being able to express a volume integral J Q q a dV in terms of boundary 
integrals is that the unit vector lie in KowAi. This is equivalent to the 
consistency of (35). Invertibility of M. is sufficient but not necessary. It is an 
open question whether, by including a growing set of vector functions {p/?} 
(Table 2 and its continuation), and necessarily the resulting growing set of 
scalars {q a }, that the row space of the resulting growing matrix Ai can be 
made to contain the unit vector b*") for any desired a. 

Note that the methods of this appendix are probably not simply general- 
izable to Riemannian manifolds with boundary. The Euclidean boundary for- 
mula (11) is known to break down in a general manifold, and unit-normalized 
eigenfunctions with exponentially-small values and gradients everywhere on 
the boundary become possible [13]. This is associated with the existence 
of trapped geodesies. However for the constant-curvature case, it would be 
interesting to search for a generalization. 
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